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Abstract 

The ADI iteration is closely related to the rational Krylov projection methods for con- 
structing low rank approximations to the solution of Sylvester equation. In this paper we 
show that the ADI and rational Krylov approximations are in fact equivalent when a special 
choice of shifts are employed in both methods. We will call these shifts pseudo %2-optimal 
shifts. These shifts are also optimal in the sense that for the Lyapunov equation, they yield 
a residual which is orthogonal to the rational Krylov projection subspace. Via several exam- 
ples, we show that the pseudo %2-optimal shifts consistently yield nearly optimal low rank 
approximations to the solutions of the Lyapunov equations. 

1 Introduction 

Let A G R nxn , B G R mxm and Y G R nxm be given matrices. Then, the Sylvester equation 
for the unknown matrix X G R nxm is given by 

AX + XB + Y = 0. (1) 

The equation (1) has a unique solution if and only if Xi(A) + Xj(B) / for i = 1, . . . , n 
and j = 1, . . . , m. A special case of the Sylvester equation is the Lyapunov equation, where 
B = A* and Y = Y* > 0. Both the Sylvester and Lyapunov equations are an important 
tool in the analysis of asymptotically stable linear dynamical systems of the form 

x(t) = Ax(t) + bu(t), y(t) = c*x(t), (2) 

where A G R nxn and b, c* G R n . In (2), x(t) G R n , u(t) G R, y(t) G R, are, respectively the 
state, input, and output, of the underlying system. While the cross gramian X of (2) solves 
the Sylvester equation 

AX + XA + be* = 0, (3) 
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the controllability gramian "P and the observability gramian Q solve the Lyapunov equations 

AT + V A* + bb* = and A*Q + QA + c*c = 0, (4) 

respectively. These three gramians are of fundamental importance especially in the concept 
of model reduction, see [1]. In what follows, we will mainly focus on the Sylvester equation 
(1) where Y is rank-1; hence our discussion already contains the Lyapunov equations as a 
special case. 

The standard direct method for solving (1) is due to Bartels and Stewart [3]. However, this 
method requires dense matrix operations such as the Schur decomposition; thus is not appli- 
cable in large-scale settings. For large-scale settings, iterative methods have been developed 
that take advantage of the sparsity and the low-rank structure of Y. The two most common 
ones are the Alternating Direction Implicit (ADI) Method ([30, 8, 9, 27, 17, 37, 33, 21, 31, 
43, 39, 44]) and the (rational) Krylov projection methods ([22, 25, 13, 24, 32, 38, 14, 2]). 

The ADI method was first introduced by Peaceman and Rachford [29] for solving parabolic 
and elliptic PDEs, and was later adapted to solving the Sylvester equation by Wachspress in 
[43]. It is a fixed point iteration scheme for approximating X . Given two sequences of shifts 
{ax, ct2, • • • , a r , . . . }, {/3i, 02, • • • , Pri • • • } C C and an initial guess Xq, the ADI iteration for 
(1) proceeds as follows : 

Xi =(A - oil) (A + fclY^Xi^B - PJ){B + aj)- 1 (5) 
-(ai + PiXA + fcl^YiB + ail)- 1 . (6) 

The performance of the ADI iteration depends heavily on the choice of shifts used in 
the iteration. Several schemes have been developed for making asymptotically optimal shift 
selections if some information is known about the boundaries of the numerical range of A, 
and B. See [42, 43, 35, 36, 31] and the references therein for further details on the shift 
selection problem in the ADI iteration. 

A closely related method to the ADI iteration is the rational Krylov projection method 
(RKPM). In the RKPM, the Sylvester equation AX + XB + Y = is projected onto 
the rational Krylov subspaces /C£ at (A, b, a) = span{(<7i/ — A)~ 1 b, . . . , (ay J — A) _1 b} and 
K7 r at {B* ,c, p) where cr= {a\, . . . ay}, and p, = {fii, . . . , fi r } are the sets of shifts used to 
construct the respective rational Krylov spaces and v denotes the conjugate of v. See [6] 
for further details regarding ZC£ at (A, b, a), and constructing an orthonormal basis via the 
rational Arnoldi iteration. Let Q r and U r denote the orthonormal basis for /C£ at (A, b, a) 
and fC^ (B* , c, p) . Then, the RKPM approximation is constructed by first solving 

Q* r AQ r X r + X r U*B*U r + Q* r bc*U r = (7) 

and then approximating X by Q r X r U* . The solution of the projected Sylvester equation 
(7) is very cheap. Like the ADI method, the RKPM method also relies heavily on a good 
choice of shifts to produce accurate results. In the next section we will derive results that 
show for a certain choice of shifts, the RKPM and ADI methods are indeed equivalent. 

Since in almost all applications, the quantities A, B, b, and c are real, we will assume 
that the set of shifts a and fi are closed under conjugation so that the approximants are 
real as well. This will guarantee that the orthonormal bases Q r for ZC£ at (A, b, a) and U r for 
fC T r at (B*, c, p) can be computed to be real as well. 
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2 Equivalence of the ADI and Rational Krylov Pro- 
jection Methods for pseudo-%2 optimal points 

In this section, we present our main results illustrating the connection between the ADI and 
RKPM. Since the discussion requires the concept of %2-optimal points for model reduction, 
we first briefly review the %2 approximation problem. 



2.1 Optimal %2 model reduction 

For a full-order model as given in (2), the model reduction problem seeks to construct a 
dynamical system 

x r (t) = A r x r (t) + b r u(t), y r {t) = c* r x r (t) (8) 

of much smaller dimension r <C n, with A r 6 W XT and b r , c* 6 M r such that y r (t) approxi- 
mates y(t) well for a wide range of inputs u(t). The reduced-model in (2) is usually obtained 
via state-space projection: Two matrices V r , W r £ ]R nXT ' are constructed with W*V r = I r 
to produce 

A r = W?AV r , b r = W?b, and c r = V r c (9) 

One can measure the quality of the approximation using the concept of transfer function. 
By taking the Laplace transforms of (2) and (8), one obtains the transfer functions H(s) = 
c(sl — A) _1 b and H r (s) = c r (sl r — A r ) _1 b r , respectively. Hence, one can consider model 
reduction in terms of these transfer functions as approximating a degree-n rational function 
H(s) with a degree-r one H r (s). For more details on model reduction of linear dynamical 
systems, see [1]. 

In this paper, we focus on the %2- n orm to measure accuracy of the reduced-model. The 
H2 optimal model reduction problem seeks to construct a reduced system as in (8), so that 
H r (s) minimizes the TI2 error over all stable linear dynamical systems of the form (8), i.e. 



\H-H r \\ H2 = min \\H - H r \\ H2 (10) 

deg(H r )=r 



where 



\H-H r \\ n2 = J \H(iu) - H r (ioj)\ 2 du?j 
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Several methods have been introduced to solve (10); see, for example, [34, 23, 45, 46, 28, 
19, 18, 41, 11, 16, 4, 5], and the references therein. Since the optimization problem (10) is 
nonconvex, the common approach involves finding reduced-order models satisfying the first- 
order necessary conditions of T^-optimality. The next theorem states the interpolation-based 
necessary conditions for Hi optimality introduced by Meier and Luenberger [28]. 



Theorem 1. [28, 19] Given a full-order system H(s) of order n, if H r (s) = ^ ^3^- is an 
%2 optimal approximation to H(s), then 



s—Xi 
i=l 



H(—Xi) = H r (—Xi) for i = l,...,r, and (11) 
H'(-Xi) = H' r (-Xi) for i = l,...,r (12) 
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A reduced-order model which satisfies the %2-optimality conditions can be obtained by 
using the Iterative Rational Krylov Algorithm (IRKA) of Gugercin et. al. in [19]. However, 
in this paper we will focus on satisfying only (11) (without the derivative condition). We will 
call these interpolation points pseudo %2- optimal points to emphasize that they only satisfy 
a subset of the optimality conditions. In terms of the projection framework (9) for model 
reduction, this corresponds to finding interpolation points a = {cii, . . . ,a r } and choosing, 
in (9), V r = W r = Q r where Q r is an orthonormal basis for the rational Krylov subspace 
/C£ at (A, b, a) so that {Ai, . . . , A r }, i.e. the eigenvalues of A r = Q^AQ r , become the mirror 
images of the interpolations points cr = {o~i, . . . , o~ r }, i.e. 

A(A r ) = \(QjAQ r ) = -a. (13) 

The emphpseudo %2-optimal points interpolation points can be computed iteratively in a 
manner similar to IRKA [19] as done in [20] for port-Hamiltonian systems. 

2.2 The ADI Iteration and Rational Krylov Projection Method 

The main theorem requires the following lemma, which connects the ADI approximation for 
the Sylvester equation with rational Krylov subspaces. This extends an earlier result by 
Li and White [26] which establishes a similar connection for the the case of the Lyapunov 
equation. 

Lemma 1. Let Y = be* , where b 6 1" and c G M. m . Let {ai, . . . ,<7 r } and {/xi, . . . , /j, r } 
be two collections of shifts that satisfy 5ft(/Uj), ;K(o"j) > for i = l,...,r. Suppose X r is 
the approximate solution to the Sylvester equation (1) obtained by applying the pair of shifts 
ai = —at and /3j = — m in the ADI iteration (5) for i = 1, . . . ,r with Xq = 0. Then there 
exist L r £ C nxr and R r £ C mxr such that X r = L r R* and colspan(L r ) C /C™*(A, b, fx) and 
colspan(Rr) C K™*(B* ,c,&) 

Proof. The proof is given by induction on i, the iteration step. First note that for i = 1, X\ = 
(/xi + 01) (A- mI)- l bc*(B-oiiy l , so let L x = [(mi + <ri)(A - ^xl)~ l b] and Ri = [(B* - 
cfi/) _1 c]. Then L\ and R\ clearly satisfy the hypothesis and X\ = L\R\. Now suppose that 
the statement holds for Xi. Then, for j = 1, . . . ,i, the j th column of Li is T^\A)b, where 
Tp\\) is a proper rational function that lies in the span of { A J L ^ 1 , • • • , tz^t}- Similarly, the 

j' th column of Ri is S^(B*)c, where (A) lies in the span of { ^j^r , . . . , }. Therefore 
Xi + \ can be written as 

X l+1 =(A + a l+1 I)(A - ^ i+l I)- l LR*{B + ^ i+1 I)[B - cj i+l B)- 1 
+ (/i i+ i + a i+ i){A - in+iI)- l bc*{B - (J i+1 I)- 1 

i 

= ^2(A + a i+l I){A - /^ +1 /)- 1 r i ' ) (A)6c*^ ) ( J B)( J B + fi l+1 I)(B - a^I)- 1 
3=1 

+ + a l+ i){A - W+ iI)- 1 bc*( J B - a i+1 I)- 1 

For j = 1, . . .i, let the j th column of L i+ i be (A + a i+1 I)(A - / u i+ i/)" 1 T. ( '' ) (A)b and let the 
(z + l) th column be T^ 1 ^(A)6 = (fi i+1 + o- i+ i)(A- fi i+ il)~ 1 b. Then clearly colspan(Lj + i) C 
/C^\(A,6,/i). Similarly, let (B* - ai +1 I)- x (B* + fn+i^S^ (B*)c be the j th column of 
fli+i for j = 1, . . . ,i, and S^\b*)c = (B* - a i+1 I)- l c be the (i + l) th column. Then 
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colspan(i?j + i) C /C^\(B*, c, a). Finally, we note that by construction, Xi + i = Li + iR* l+l . 

□ 

Next, we give our first main result showing that the approximate solution of the Sylvester 
equation (1) by ADI and RKPM are indeed equivalent when the shifts are chosen as pseudo- 
V.2 optimal points. This result applied to the special case of Lyapunov equation was first 
presented at the 2010 SIAM Annual Meeting [15] then later published independently in [13]. 
Our new result here, on the other hand, is more general than both [15] and [13] since it 
tackles the case of Sylvester equation and includes the Lyapunov equation as a special case. 
Moreover, while the proof given in [13] for the special case of Lyapunov equation makes use 
of a novel connection between the ADI iteration and the so-called Skeleton approximation 
framework first developed in the work of Tyrtyshnikov [40] , the proof we provide here for the 
more general Sylvester equation case is given directly in terms of rational Krylov interpolation 
conditions, and in that sense is simpler. 

Theorem 2. Given the Sylvester equation (1) with Y = be* , where b E R n and c E M m , 
let Q r E M nxr be an orthonormal basis for the rational Krylov subspace JC™ {A, b, <r) and 
let U r E M mxr be an orthonormal basis for the rational Krylov subspace K^{B*, c, a) for 
a set of shifts cr = {a\, . . . , o~ r } where 3?(<Tj) > for i = 1, . . . , r. Let X r E W xr solve the 
projected Sylvester equation 

Q* r AQ r X r + X r U*BU r + Q*bc*U r = 0, (14) 

and let X r E ]R nxm be computed by applying the shifts on = —o~i and Pi = —o~i to exactly r 
steps of the ADI iteration (5) for i = 1, . . . ,r. Then X r = Q r X r U* if and only if either 
\{Q* T AQ r ) = -{(Ti, ...,a r } or \{U*BU r ) = -{a u ...,a r }. 

Proof. (<!=) First suppose that \{Q*AQ r ) = — {<7i, . . . ,a r }. The proof remains the same if 
we instead suppose that X(U*BU r ) = — {ai, ■ ■ ■ , oy}. Let A r = Q*AQ r , and b r = Q*.b, 
B r = U*BU r , and c r = U*c. Note that after we apply r steps of the ADI iteration 
with the set of shifts at = (3% = — o~i to the projected Sylvester equation (14), we obtain 
the exact solution X r . since X(A r ) = —{a±, . . . ,a r }. By Lemma 1, at the rth step of the 
ADI iteration X r = L r R* where L r = [T^(A r )b r , . . . , T^\A r )b r ) where T^(A r )b r are 
rational functions that lie in ZC£ at (A r , b r , a). Similarly R r = [S ( - 1 \B*)c*,...,S^\B*)c*] 
where the S^(B*)cr are rational functions that lie in ZC£ at (-B*, c*, cr). Furthermore, for the 
same shifts, on = Pi = —o~i for i = 1, . . . , r, applied to r steps of the ADI iteration on the 
full Sylvester equation (1), we have X r = L r R* and L r = [T^(A)b, . . . ,T^(A)b] and 
Rr = [S^(B*)c,...,S^(B*)c]. Thus it is sufficient to show that Q r L r = L r and that 
U r R r = R r . Without loss of generality consider just the former equation. This, in turn, 
amounts to showing that Q r T®(A r )b r = T®(A)b. IiTi(A)b are a set of orthogonal rational 
functions that span /C£ at (A, b, a), then it is sufficient to show that 

Q r Ti(A r )b r = Ti(A)b. (15) 

Equality (15) follows readily from the interpolation properties of the Galerkin projection, 
which we show below. First, note that due to the interpolation properties of the Galerkin 
projection, Q r [pil r — A r )~ x b r = (crjJ — A)~ x b. Let V r = [{ail - A) _1 b. . . {a r I - A) _1 b]. 
Then, for some x E M r , 

V r x = Ti{A)b = Qr[{ail r - A r )~ l b r . . . {a r I r - A^b^x = Q r Ti{A r )b r , (16) 
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which proves (15). 

(=>) Let X r be the solution of 

Q* r AQ r X r + X r U*BU r + Q* r bc*U r = (17) 

where Q r is an orthonormal basis for /C£ at (A, b, a) and U r is an orthonormal basis for 
lC r r at (B*, c, a). Suppose that Q r X r U* = X r . Let X r be the approximate solution of (17) 
resulting from applying the shifts oti = fi% = —Oi for i = 1, . . . , r to exactly r steps of the ADI 
iteration (5). By the interpolation result given in the proof above, Q r X r U* = X r . It follows 
from the assumptions that, Q r X r U* = Q r X r U*, so X r = X r . But this means that X r 
solves (17), and so either X(Q*AQ r ) = —{ax, . . . , oy} or X(U*BU r ) = —{ax, ... , oy}. □ 

Remark 1. This theorem shows that the ADI approximation for the Sylvester equation is 
equivalent to lifting the solution of the projected Sylvester equation back to the original di- 
mension when either X(Q*AQ r ) = — {ax, ■ ■ ■ , oy} or X(U*BU r ) = — {ax, ■ ■ ■ , ay}; hence the 
two most common approximation methods for solving a Sylvester equation is indeed equiva- 
lent for these special shift selection. Recalling pseudeo-%2 optimality condition (13), for a 
given r, these special shifts are indeed exactly the pseudo-%2 optimal shifts for a dynamical 
system Hx(s) = Zx(sl — A)~ l b or i?2(s) = ^(sJ — B*)~ l c* where z\ and Z2 are vectors of 
appropriate sizes. 



2.3 Orthogonality in the case of Lyapunov equation 

The parameters for which the ADI iteration and the rational Krylov projections coincide also 
satisfy orthogonality conditions on the residual for the special case of the Lyapunov equation 

AX + XA* + bb* = (18) 

For a given approximation X r to the solution X , define the residual R as 

R = AX r + X r A* + bb* . (19) 

The following result was first given in [13]. Here we present a new and more concise proof 
of the orthogonality result in terms of the special interpolation properties of the pseudo 
%2-optimal shifts. 

Theorem 3. Given AX + XA* + bb* = 0, let X r G W xr solve the projected Lyapunov 
equation 

Q*AQ r X r + X r Q*AQ r + Q* r bb*Q r = 0, 

where Q r is an orthonormal basis for the JC™ (A, b, cr) with a = {a±, . . . ,a r } Let X r = 
Q r X r Q*. Then Q*R = if and only if X(Q*.AQ r ) = —{ax, • • • , oy} where R is the residual 
defined in (19). 

Proof. (=>) Suppose that Q*R = 0. Multiplying (19) with Q* from the left and then 
transposing the resulting equation leads to 

AQ r X r + Q r X r Q* r A*Q r + bb*Q r = 0. (20) 

Let A r = Q*AQ r = TAT -1 be the eigenvalue decomposition of A r where A = diag(Ai, . . . , A r ). 
Plug these expressions into (20), and right multiply by T~* to obtain 

Q r X r T~*A* + AQ r X r T~* + bb*Q r T* = (21) 
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Let Q be the v entry of b*Q r T * . Then it is straightforward to show that the v column of 
Q r X r T* must be (-A<7 - A)~ l bd. Thus, it follows that K,f\A,b,<T) = K,f\A,b, -A), 
where A = {Ai,...,A r }. Since both sets <x and A are closed under conjugation, after an 
appropriate reordering, we obtain Oi = — A$. 
(<=) Observe that 

A r X r + X r A; + Q* r bb*Q r = =>■ (22) 
A r X r T _ * + X r T~*A* + Q* r bb*Q r T * = 0. (23) 

Thus, the it/i column of X r T~* is (-A»/ r - A r ) -1 Q*6<i. But since Q r is an orthonormal 
basis for /C£ at (A, b, a), and Aj = — <7j, this means 

Qr{-\il r - A r )~ l Q* r bQ = {-Xil - Ay l bQ = (Q r X r T-*)e u (24) 

where e, is the it/i unit vector. Thus, 

Q r X r T~*A* + AQ r X T T * + bb*Q r T * = 0, (25) 

which implies 

Q r X r Q*A*Q r + AQ r X r + bb*Q r = 0. (26) 
Transpose this last expression and use the fact that Q*Q r = Ir to obtain 

Q*Q r X r Q*A* + Q* r AQ r X r Q* r + Q* r bb* = Q*R = 0, (27) 
which is the desired result. □ 

Remark 2. As we have previously noted, in almost every practical situation, one would 
choose a set of shifts cr which is closed under conjugation. But even for the cases where this 
assumption on a does not hold, Theorem 2 holds as is, and Theorem 3 applies with a slight 
modification. To wit, Aj = U{ for i = 1, . . . , r. 



3 A numerical study on using the pseudo-%2 opti- 
mal points as the ADI shifts 

Having shown that using the pseudo-%2 optimal points in the ADI iteration for the Sylvester 
equation is equivalent to applying RKPM and that the pseudo-%2 points leads to an orthogo- 
nality condition in the case of Lyapunov equation, the natural question to ask is what quality 
of approximation the pseudo-%2 optimal points have as ADI shifts. We will briefly inves- 
tigate this issue in this section. However, we emphasize that the purpose of our numerical 
results is not to advocate employing the pseudo %2-optimal shifts in the ADI iteration or in 
the RKPM. This would be a costly numerical method for approximating Sylvester equations 
since obtaining the pseudo %2-optimal shifts already requires solving several linear systems. 
Our numerical results are meant to illustrate the unique quality of these shifts compared 
with other choices of shifts that do not share the ADI-RKPM equivalency property. 

We used three benchmark models in our numerical simulations: The CD Player model 
with n = 120, the EADY model with n = 598, and the Rail Model with n = 1357. The first 
two models are described in detail in [12] and the Rail model in [7]. For all three models, 
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we compute a rank r approximation to the solution of the the Lyapunov equation (18). The 
exact and approximate solutions are denoted by X and X r , respectively. The Rail model 
has multiple inputs; thus for this model we only use the first column of the input matrix. 
We use three different approximation methods for each model: 

• Method 1: The RKPM is applied to the a sequence of shifts that alternates between 
and oo. The resulting subspace is generally referred to as the extended Krylov subspace. 
Its application to RKPM was first introduced by Simoncini in [32]. 

• Method 2: The RKPM is applied using r pseudo-%2 optimal shifts; or equivalently r 
steps the ADI iteration is applied using r pseudo-%2 optimal shifts. 

• Method 3: The r-steps of ADI iteration are applied where the ADI shifts are chosen 
via Penzl's heuristic method [30]. 

The quality of the resulting approximations from each method are compared using the 

lljf X II2 

relative error in the 2-norm, i.e. r — . Figure 1 shows the relative errors for the EADY 

model as r varies from 1 to 50 together with the minimum possible error, i.e. where 7Tj is 
the i th singular value of the true solution X. Note that for a given r, the pseudo-%2 optimal 
shifts perform remarkably well, almost matching the best low-rank approximation given by 
the singular value decomposition. For a selected number of r values, these numbers are also 
tabulated in Table 1 further illustrating the effectiveness of the pseudo-%2 points as ADI or 
RKPM shifts. Similar results for the CD player model are shown in Figure 2 and in Table 2 
and for the Rail Model in Figure 3 and Table 3 illustrating that the pseudo %2-optimal shifts 
produce a nearly optimal rank r approximation in several cases. Indeed, this phenomenon 
was recently explained by Breiten and Benner in [10] , where they show that the %2 optimal 
shifts are optimal in a special energy norm related to the Lyapunov equation; for details we 
refer the reader to [10]. 




5 10 15 20 25 30 35 40 45 50 
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Figure 1: Relative error as r varies for the EADY Model 



8 



Table 1: Comparison of Methods with the EADY model order 598 

\\X — X r \\2 



X 2 



r 


7I>+1 

2.31 x 10" 




Method 1 


Method 2 


Method 3 


10 


-4 


8.42 x lO" 2 


1.28 x 10~ 3 


1.96 x lO" 1 


20 


3.38 x 10" 


-7 


5.73 x 10- 2 


4.99 x 10~ 7 


1.13 x lO" 1 


30 


7.63 x 10" 


-9 


1.09 x 10~ 3 


8.46 x 10~ 9 


8.54 x 10~ 2 


40 


4.83 x 10" 


-10 


1.93 x 10~ 4 


1.06 x 10~ 9 


9.70 x 10-3 




10 20 30 40 50 60 70 80 

Rank of Approximation 



Figure 2: Relative error as r varies for for the CD Player model 
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5 Conclusions 

In this paper we presented a new result that solidifies the connection between the ADI 
iteration and rational Krylov projection methods for solving large-scale Sylvester equation. 
We have shown that for one-sided projections, the two methods are indeed equivalent for 
a special choice of shifts called pseudo-%2 optimal shifts, so-called because they partially 
satisfy first-order necessary conditions for %2 optimal model reduction. These shifts are also 
optimal in the sense that they produce an approximation with a residual orthogonal to the 
rational Krylov projection subspace in the case of Lyapunov equation. 
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